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0> Abstract 

1 I In this paper, we present new, explicit, volume-preserving vector fields 

for polynomial divergence-free vector fields of arbitrary degree (both pos- 
itive and negative). The main idea is to decompose the divergence poly- 
i nomial by means of an appropriate basis for polynomials: the monomial 

basis. For each monomial basis function, the split fields are then identified 
^ by collecting the appropriate terms in the vector field so that each split 

a vector field is volume preserving. We show that each split field can be 
I I integrated exactly by analytical methods. Thus, the composition yields a 

volume preserving numerical method. Our numerical tests indicate that 
( the methods compare favorably to standard integrators both in the quality 

^ of the numerical solution and the computational effort. 

Keywords: Geometric integration; volume preservation; splitting 
methods 

0^ 

lO 1 Introduction 

O 

CN Diver gence-free vector fields occur naturally in incompressible fluid dynamics, 

" , I and preservation of phase-space volume is also a crucial ingredient in many, if 

^ not all, ergodic theorems. Preservation of volume by a numerical method for 

differential equations is thus a desirable property in the study of dynamical 
system. Nevertheless, designing volume preserving numerical integrators is a 
hard task, as the space of divergence-free vector fields seems to be too large. [S] 
proved that no standard numerical method, for example a Runge-Kutta scheme, 
is volume preserving for all such vector fields. Thus, the design of efficient 
methods that preserve volume is still a standing open problem in geometric 
integration [7]. 

Despite no-go theorems for volume preservation within the class of "stan- 
dard" methods (see also the recent results in [3], [T], where it is proved that no 
B-series method can be volume preserving for all possible divergence-free vector 
fields), it is known that volume preservation can be achieved, either restricting 
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the class of vector fields or using methods other than B-series. [T] have shown 
that volume-preserving B-series can be obtained if the vector field has some 
specific dependence (for example /i = /i(x2),/2 = /2 (xa ),..., /„ = /n(a;i), 
or if the variables naturally decompose into two sets, x = [y,z]-^, obeying 
y = g(z),z = h(y)). These two examples of vector fields correspond to what 
we call off- diagonal. These are generally easier to treat and other methods will 
be described in the sequel. 

Volume-preserving maps can be constructed techniques other than methods 
that possess B-series, for instance using generating functions [101 HSl US] ■ This 
technique involves evaluating definite integrals of the vector field, and, in ad- 
dition, the method does not preserve fixed points. Another technique is based 
on splitting methods. One of the earliest volume-preserving splitting methods 
is indeed the splitting method by [4j, decomposing the vector field into the 
sum of essentially 2-dimensional Hamiltonian fields, which are then solved by a 
(typically implicit) symplectic method. 

Because of the difficulty of addressing the general space of divergence-free 
vector fields, recent efforts have concentrated to smaller, yet still interesting, 
functions spaces, for instance the space of polynomial fields. An earlier paper on 
splitting polynomial vector fields is by [6] . That paper had some discussion of the 
divergence-free case, but mainly dealt with the Hamiltonian case. Investigations 
of the Hamiltonian case, which involves expressing a scalar polynomial of degree 
d in n variables as a sum of functions of fewer variables, have shown that good 
splitting methods exist, but that finding and analyzing them (especially for 
general n and d) is very difficult [21 [6l E] . The volume-preserving case, which 
involves n polynomials subject to the divergence-free condition, is even harder, 
although there is a conjecture by [B] that they can be expressed as a sum of 
n + d shears, each a function of n — 1 variables. The case of linear and quadratic 
divergence-free vector fields was studied in detail in [Q, where several explicit 
volume-preserving splitting methods were introduced. In that paper, two main 
classes of methods were considered: a) methods that distinguish the diagonal 
and off-diagonal part and b) methods that do not. By diagonal part we mean 
all the terms of the vector field such that Xi depends on a;^, for i — 1, . . . ,n. 
Similarly, the off-diagonal part refers to all the terms of the vector field such 
that Xi does not depend on Xi, i = l,...,n. Furthermore, for the class a), 
several explicit schemes that treated the diagonal part and the off-diagonal part 
separately were introduced and tested numerically. Numerical tests indicated 
that methods that treated the diagonal part by splitting it in terms treated by 
exponentials had a smaller error than methods splitting in shears. 

In this paper we present a new approach that allows us to develop ex- 
plicit volume-preserving methods for arbitrary polynomial divergence-free vec- 
tor fields, including those with negative degree. The first main insight is to 
expand the divergence equation, rather than the vector field, in the monomial 
basis. For each monomial basis element, we identify the elements in the vector 
field associated to it to construct a divergence-free elementary vector field. The 
second insight is to recognize that the elementary divergence-free vector fields 
can all be treated by the same formalism and therefore can be solved explicitly 
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by elementary analytical methods. The split fields are then composed to obtain 
explicit first order method and second order method (by symmetrization) . The 
resulting composition method is thus explicit and volume-preserving. Being ex- 
plicit, the proposed method are computationally efficient and possess excellent 
qualitative properties. We believe that thisis a consequence of volume preser- 
vation solely, as the methods are not necessarily time-reversible nor self-adjoint 
(for instance, the first order method is neither). 

The paper is organized as follows. In Section [2] we review some background 
and introduce notation. In Section [3] we present the monomial basis for polyno- 
mial volume-preserving vector fields and prove that the vector fields associated 
to the basis elements can be integrated exactly. The case of polynomials of 
negative degree is treated in Section [5j In Section |6] we test the application 
to volume-preserving cubic Stokes flows. Finally, Section [7] is devoted to some 
concluding remarks and some directions for future research. 

2 Background and notation 

We consider the ordinary differential equation 

x-f(x), x(0)=xo, (1) 

where x e M" and f : M" ^ M", f(x) = [/i(x), . . . ,/„(x)]^, is subject to the 
divergence-free condition 

n 

V-f = ^9,J,(x)=0. (2) 

i=l 

An arbitrary vector field f(x) can always be decomposed as 

f(x) =f'^'''S(x)+f°ffdiag(^) 

where, component-wise, ff^^^{x) is the collection of terms in /i(x) that depend 
on Xi ( i.e. 9a;;/f"^^(x) ^ 0). Similarly, P^'^'^^s^x), is given, component-wise, by 
the collection of terms in fi that do not depend on Xi, i.e. dxif°^'^^^^{^) = 0. 
We refer to f^'^^s and f°ffdiag^^-j ^^^^ diagonal part and the off-diagonal part 
of the vector field f , respectively. 

From the definition of divergence, it is clear that only the coefficients of the 
diagonal part are directly involved in the divergence- free condition ([2| , therefore, 
vector fields with zero diagonal part are automatically divergence-free. 

The off-diagonal part of a vector field is generally easier to treat by volume- 
preserving methods. The method of splitting in canonical n-shears is generic 
(for each i = 1. ..n, we solve for Xi, while Xk = 0, k i). For polynomial 
systems, it is possible to construct splittings in lower-triangular systems (for 
each j = 1 . . . n, we solve for Xi depending only on . . . , Xi-i, under a suitable 
permutation of the indices), see [5]. 
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Thus, to obtain volume-preserving methods, special care has to be taken 
when splitting the diagonal part, because splitting across the conditions for 
zero-divergence will yield methods that are not volume-preserving. The goal 
is to split the diagonal part in divergence-free vector fields that can be solved 
exactly. These observations motivate the following assumption. 

Assumption 1 Unless otherwise stated, we will assume that the given vector 
field is diagonal only (i.e. f = f'^^^s) and that the off-diagonal part is zero 

^foffdiag _ Qj 

In this paper we consider the case when the functions /i(x) are polynomials. 
To treat the general case with n variables and given degree d, it is convenient 
to introduce a multi-index notation. 

Definition 1 Letj = {ji, ■ . ■ , jn) G f^" be a multi-index and let = x-f^x-!^ ■ ■ ■ xii 
(monomial) . In addition, denote by\]\ = Ji-f • ■ •-\-jn the degree of the monomial. 

Earlier studies of splitting methods for polynomial fields have focussed on 
homogeneous polynomials. Assume that the /i(x) are homogeneous polynomials 
of degree d in the variables xi,X2, . . . ,Xm i-C. 

ii = /j(x) = ^ aj xj, i = l,...,n. (3) 

\i\=d 

The case when /i(x) is not homogeneous is easily treated by a further splitting 
in homogeneous terms. Define N{n,d) = In [B] it was shown that 

N{n, d) is the number of coefficients aj = a*^ (for a fixed value of i) of 
an homogeneous polynomial fi of degree d in n variables. The divergence-free 
condition is a homogeneous polynomial of degree d—1 because of derivation. All 
the coefficients of this polynomial must be identically equal to zero, because of 
the arbitrariness of the variables. This gives N{n,d— 1) conditions for volume 
preservation. In the polynomial case, there are N{n, d—1) divergence-free 
conditions on the coefficients. In order to be volume preserving, a splitting 
must satisfy one or more conditions on the coefficients. Our approach is inspired 
in part from that argument. The new idea is to consider a polynomial vector 
field and to look at a monomial-basis expansion of the divergence function. 
Each basis element will correspond exactly is associated to a condition on the 
coefficients of the vector field. As long as we split according to basis terms, we 
are guaranteed to obtain "elementary" monomial fields that are also divergence 
free. We show that all these "elementary divergence-free vector fields" possess 
n — 1 integrals in evolution and are therefore integrable. Furthermore, we give 
the explicit analytical solution. This allows us to generate splitting methods for 
arbitrary divergence-free polynomial fields. 

We illustrate the general idea by discussing a simple 2-dimensional example. 
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Example 1 Consider the vector field 

xi = aj xl + x\x2 + a\ x\x% 

2 2 I 2 2,23 
X2 = (Xj X-^X2 + Clk X1X2 + Ct] X2 

The divergence of this vector field is 

p{x-i_,X2) = (Saj + a?)xl + {2al^ + 2a^)x-i_X2 + (ai + M\)x\. 

The divergence free condition becomes the set of equations 

Saj + a| = 0, Ok + Ok) "^i + 3a? = 0. 

As argued above, for each equation we obtain a divergence-free split vector fields 

Xi = aj x\ Xi = a]^ x\x2 xi = aj X1X2 

±2 = a? xlx2 ±2 = a| xixl, ±2 = a? x^. 

Surely, each of these vector fields can be integrated. For instance, in the first 
vector field, one can solve for x\ then substitute into the second equation to 
solve for X2 ■ In the third vector field, the procedure is similar, just interchange 
the role of xi and X2- The second vector field can also be integrated, though 
not exactly by the same method. The situation becomes far more complicated 
for several variables and higher order polynomial vector fields, so that, at first 
glance, this procedure does not yield a method that can be easily generalized. 
However, observe that the three vector fields can be written as 

Xi — aj xi ■ xl xi — al.xi{xiX2) Xi = al x\ ■ x\ 

±2 = (^X2- xl ±2 = Ok X2{XlX2) ±2 = X2X%, 

i.e. each vector field obeys an equation of the type Xi = CiXi{x{^ x^2) where 
x{^x^ are the first, the second and third monomial in p{xi,X2). In particular, 
if {x{^ x-2){t) is known, a;, can be obtained by integration. 

In what follows, wc shall sec that this argument is generic: it yields for any 
number of variables and any degree of the vector field. We shall therefore 
identify monomials appearing in the divergence-free condition and solve for them 
explicitly. When the monomials, as function of time, arc known, the remaining 
variables follow by integration of linear equations with variable coefficients. 

3 Elementary divergence-free vector fields: the 
polynomial case 

Consider a polynomial vector field 

d 

Xi = /i(x) = alx^, 
|k|=i 



i = l,...,n, 



(4) 
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of degree d. We assume that this polynomial field is diagonal in the sense 
described earlier. Diagonality implies some (mild) restrictions on the nonzero 
coefficients of namely that, for a given index k, one has aj^, ^ provided 
that k, ^ 0. 

Rather than looking at the vector fields, it is useful to focus on the divergence 
polynomial and reconstruct the vector field from the divergence. 

Lemma 1 Consider the (diagonal) polynomial vector field Q. Assume that 
the vector field is divergence-free. Let 

|j|=o 

be the divergence o/Q. For each multi-index], consider the polynomial field 
ii = ftj+e,^'^''' 1 < i < n. (5) 

where the canonical unit vector in M" with 1 in the ith position and else- 
where. The polynomial field ^ is divergence free. Moreover, Q can be uniquely 
split in sum of vector fields of the form ([s]) . 

Proof. Since monomials are basis for the set of polynomials of degree d, 
zero-divergence of Q implies that 

Pi=0 |j|=0,...,d-l, 

where pj is the coefficient of = ■ ■ ■ x;^" in the divergence polynomial. For 
any multi-index j, a contribution to the coefficient pj comes from dx.fi{x), for 
each 1 < i < n. If the fi{x) are polynomials, such a term originates from deriva- 
tion of Xixj = xj+''' with respect to Xi. This corresponds to the term a'-_^g.XiX^ 
in fi (x) . This procedure picks up all the terms in fi (x) contributing to and 
the condition pj = guarantees that ([5| is divergence free. The uniqueness 
comes from the fact that J x-i dxi ~ jT^xJ+'^i + c/i(x), where dxih{x) = 0, i.e. 
h is an arbitrary function of the remaining variables and does not contribute to 
the diagonal part of the vector field. ■ 



In the sequel, we will often use the short-hand notation 

x-Fj(x) (6) 

for ([5]). Each of these divergence-free vector fields is associated to a monomial 
basis element, and will be called an elementary divergence-free vector field (in 
short, EDFVF). 
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4 Properties of the EDFVF 

Without loss of generality, we can consider a EDFVF Fj of the form 

(for simplicity, we have dropped the dependence of a on the actual index j). 
The divergence free condition amounts to the algebraic relation 

a^(j + l)=0, 1 = (!,..., 1)^. 



Theorem 1 (Integrability of monomial EDFVF) The function 1 1 = xj+^ 
is an integral of the EDFVF ([7| . Moreover, if h £ M" is any other vector 
orthogonal to a, then is also an integral of In particular, it follows 

that has n — 1 independent integrals of motion, hence the EDFVF ([?]) is 
integrable. 

Proof. Let b be an arbitrary vector orthogonal to a, i.e. such that a-^b = 
(the case bi = j + 1 is just a particular choice of b). We have 

— x*' = ^ hx.Ti^^''' = ^ a.b.x.x.^x^^''' 

i=l i=l 
n 

= a,hx^+^ = (a^b)xj+'' = 0, 

1=1 

because of the orthogonality condition. Furthermore, since a e M", we can 
find further n — 2 vectors orthogonal to a and to j + 1 — bi, i.e. vectors 
b2,...,b„_i (each corresponding to the integral Ik = ^ 0), such that 
{a, bi, b2, . . . , b„_i} is an orthogonal basis of M". These integrals are function- 
ally independent, i.e. any linear combination 

ciV/i + ...+c„_iV/„_i =0 (8) 

has only the trivial solution 

Ci = C2 = . . . = c„_i = 0. 

To show this, note that 





" (bfe)ixb--''i ■ 




(bfc)i/xi 












_ (bfc)„xb'^-''" _ 







Because of the arbitrariness of condition (|8]) is equivalent to the 

linear system 

[bi, . . . ,b„ i]c = 



7 



c = (ci, . . . , c„_i)-^, which admits only the trivial solution as the matrix columns 
are orthogonal, hence linearly independent. In conclusion, each of the vector 
bfc, fc = l,...,n— 1, will give rise to an independent integral of motion, hence 
integrability of the vector field. ■ 




Figure 1: A 3-dimensional system with j = (1, 1, 1)^, a = (— |, |, |)"^ and b2 = 
a X (j + 1) (cross product). The solution of the system lies on the intersection 
of the surfaces Ii{t) = x.^^^ (lighter surface) and l2{t) — x'^^ (darker surface). 



Figure [T] shows two such integral surfaces for a 3-dimensional system ([7| 
with a = (— |, |, i)-^, j = (1, 1, 1)-^ and initial condition (1, 1, 1)-^. The solution 
(blue thick line) lies at the intersection of the two integrals of motion, Ii (yellow, 
lighter color) and I2 (purple, darker color), the latter corresponding to b2 = 
a X (j + 1). 

Among the n—1 integrals, the first one Ii — xJ^-*- is more fundamental than 
the others, because it is a direct consequence of the divergence-free condition. 
In particular, it can be always used to transform the EDFVF ^ (which has 
nonzero diagonal part) into a new system with zero diagonal part. 

Corollary 1 The EDFVF Q is equivalent to the off- diagonal divergence free 
system 

ii = ai/i(to)x"-^+*'', i = l,...,n, 
where = s.^^^ is an integral of the system. 
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The above system is obviously ofF diagonal since the ith component of — 1 + 
is zero. 

We have demonstrated integrability of the polynomial EDFVF by providing 
the existence of n — 1 independent integral of motions. However, it is well known 
that the explicit knowledge of the integrals of motions does not necessarily give 
explicit informations about the solution. Below, we give the explicit solution 
for the EDFVF ([7|. This is based on two fundamental observations: 

• Xi{t) can be reconstructed by integration provided that x'{t) is known, 

• remarkably, the basis functions x-^ always obey the same type of solvable 
differential equation, independently of the multi-index j. 

Theorem 2 (Analytic solutions of monomial EDFVF) Consider the ED- 
FVF (j?]) associated to the monomial basis element x-^. The associate basis func- 
tion xJ(t) obeys the differential equation 



dt 



X-' = CjX^-', 



where cj is the constant 



For sufficiently small At < 



Cj = a J. 



|cjxJ(io)| 



the solution of {7) is 



for i = 1, 




' CjxHto)it - to) 
a,xj(to)(t-to) 



for Cj^O 
for Cj = 0, 



(9) 
(10) 



(11) 



(12) 



Proof. By direct computation, we have: 
d 



dt 



-xJ = 



JlXiX- 



• • • -|- jfiOjjiXjiX'^X'^ 



.2j 



The second passage follows from substituting the components of the vector 
field Fj in place of Xi and the last passage follows from the definition of cj. 
Remarkably, this differential equation is exactly the same as the first one in the 
quadratic case ^5 , therefore it can be solved in the same manner, for instance, 
by separation of variables, yielding 



xJ(0 = 



xj(to) 



l-Cjxj(^o)(t-to) 



= xJ(io) l-CjxJ(io)(t-to 



(13) 
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in the interval [to, to + At]. 

Thereafter, each of the equations ^ is of the form Xi — h{t)xi{t), where 
h{t) — aixj(t) is a known function. As for the quadratic case (see [5]), the 
equation is hnear in x^, with variable coefficients, and has solution Xi(t) = 

Xi{to)e'^*o '^^'^'^'^'^ . Hence, for cj ^ 0, 

= a;,(io)(l-Cjxj(to)(i-io)) 

where rj is defined as in ( 12 1. The case cj = follows in a similar manner. This 
completes the proof of the statement. ■ 



Example 2 Consider the divergence-free differential equation 



x\ = + ^3 sina;i C0SX2 + 4x2a;3 

±2 = — ^a;fa;2 — 52:3 COS sin a;2 — 

±3 = — ^a;fa:3 — |;X§ cosxi cosa;2 + xfa;2. 



(14) 



This system consists of a diagonal part and a off- diagonal part. In many ap- 
plications, it is common to approximate such vector fields by a polynomial field 
of low degree (like linear, quadratic, etc.). This is typically done by using a 



Taylor expansion and truncating to a certain degree d. 
the system 

X = f'^"'S(x)+f°ffdiag(^)^ 

where 



The choice d — 5 gives 
(15) 



X1X3 + X 



-2X2X3 - 

-ix2 
4-''3 



13 12 
1 " g-^l-^S 2'^1'^2'^3 

2X^X2 + j2 2^2^3 + 42;iX2X3 
" " - ~ ~ - " 2 

3 



2X^X3 + 0X^X3 + 0X2X 



foffdiagj-^-j _ 



4X2X3 
-Xi 
X^X2 



The off-diagonal part can be treated using canonical .shears. For the diag- 
onal, observe that the divergence- free condition for ( 15 ) is a polynomial in 
X3, xf , xf X3, X2X3. These correspond to the multi-indices j = (0,0,1), k = 
(3,0,0), 1 = (2,0,1), and m = (0,2,1). Each of these multi-indice s is asso- 



ciated to an EDFVF. 



obtain 



4 „2 



6 „3 



24 



and (12 I gives 



For instance, for \, one has c\ 

- — |. Assuming that the solution is known at tk, we 



xi{t) 

X2{t) 
X3{t) 



Xl{tk){l 
X2itk){l 
X3itk){l 



5 

24 
5 

24 
5 

24 



Xlitk)x3{tk){t - tk))-^/\ 

xl{tk)X3{tu){t-tk)f/\ 

xl{tk)x3{tk){t~tk)f'\ 



for t G [tk,tk+i\. A similar treatment applies to the EDFVF associated to the 
other multi-indices. 
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5 Extensions to the case of negative (Laurent 
polynomials), rational and real mult i- indices 
bases 

In this section, we investigate the extension of the theory to treat the case of 
sphtting in terms of the form x-^ with j e Z", or, more generaUy, j e M". 

We assume that the divergence of the vector field can be expressed in the 
form 

P{xi,X2,...,Xn) = ^Pix\ (16) 

where the C F" (where F = Z,M) is discrete finite set of muhi-indices. 
Lemma 2 Assume that the divergence free vector field ([T]) has components 

/,(x)=^aix', (17) 
where £ C F" (discrete set), i — 1, . . . ,n. Then, the divergence of t is a gener- 



alized polynomial of the form (16 I. Moreover, either the indices] in (16) have 
components ji ^ —1, or, if some ji — —1, then equation i does not contribute 
to this multi-index. 

Proof. Terms contributing to the divergence originate from dfi /dxi . Any x' = 



c'* • • • a;-" contributes aJ/iX^^ ' " ^'i ^ " ' ^l" = aj/^xj, to ( 16 ), j = 1 — e.j 



and it is easily seen that linear combinations of these terms yield a "general- 
ized" polynomial of the form ( [l6| . As far as the second part of the statement 
is concened, we see that ji = —1 (yielding a term of the type x~^ in the diver- 
gence) implies h = (this contribution is always identically zero). This rules 
out that a contribution to the term x~^ can come from the ith equation but it 
does not rule out such a contribution to come from another equation. ■ 



As a matter of facts, a contribution to the term from equa- 

tion i can occur if and only if a term of the type xj^ • • • log Xi ■ ■ ■ x^^ appears 
in fi, which is ruled out from the hypotesis that fi is of a Laurent polynomial 
type. 

The results in Theorem [l][2] can now be extended with minor modifications 
to treat the generalized polynomial case. 

Theorem 3 Let he given the differential equation where f component-wise 
satisfies (17). Let 

o^y^Pi^' (18) 



be the corresponding divergence-free condition. For each multi-index j £ i/, 
there exists the associated elementary divergence-free vector field 

Xi = {1 - Si^i)al^^^XiX.\ i = l,...,n, ji = -1, (19) 
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where Si^i is the Kronecker delta and I is any index such that ji = —1. The 
elementary divergence free vector fields { 19 1 can be written in the form 



(disregarding the dependance on the index j), they are volume preserving by 
construction and possess n — 1 integrals of motions, hence they are integrable 
and their solution is explicitly given as 

Proof. The proof is very similar to the monomial case. The only difference is 

the presence of the term 5i^i, taking care of excluding the contribution of the lih. 
coefficient in aj+gj (because this contributes to a off-diagonal term, see above 
lemma) for the specific j with ji = — \. ■ 



6 Numerical examples 

We present some interesting applications to volume-preserving vector fields for 
a quadratic and a cubic Stokes flow. In our literature search, we didn't find 
any relevant examples of Laurent volume-preserving fields, and we will test the 
methods on a artificial example, for the sake of illustration. 

6.1 Quadratic Stokes flow 

We consider a quadratic volume-preserving system introduced in [8j to study 
the distruction of adiabatic invariance under separatrix crossing, 

xi = —8x1X2 -I- ex3, 

X2 llx? + 3xi+x§-3, (20) 
X3 = 2x3X2 - exi. 

The system is integrable for e = 0. When e ^ Q the system is no longer 
integrable but solutions inside the unit sphere remain bounded in the sphere. 
We consider the non-integrable case and split the system into a diagonal and 
off-diagonal part. 





Xi 


= — 0X1X2 




Xi 


= ex3. 


fdiag . 




= 8X2 


foffdiag . 




= llx^-hx§-3, 




±3 


= 2X3X2 




±3 


= — exi. 



The diagonal part corresponds to the index j — (0, 1,0), which is solved by a 
single fiow as in Prop. [Sj with cj — [0, 1, 0] [— 8, 3, 2]^ = 3 (scalar product) and 
rj - i X [-8,3,2]^. 



12 



Figure 2: The quadratic Stokes flow (20) with e = 0.1. Initial condition xq = 
[0,0,0.96]^, integration time T = 500. Top left: volume -preserving method, 
order 2, implemented with stepsize h — 0.01. Top right: ode45, using step size 
control, with options RelTol=le-6 (default value le-3). Bottom left: ode45, 
using step size control, default implementation, until T — 220, the method 
becomes unstable at T = 2.2887e-|-2. Bottom right: same as bottom left, but 
letting T = 250. 



In Figure [2] we show the numerical trajectories corresponding to e — 0.1 
and initial condition Xg = [0,0,0.96]"^, for T = 500. The second-order volume- 
preserving method (top left) is implemented as described in Algorithm 1, with 
a constant step size h = 0.01, and gives a nice and bounded trajectory. The 
first-order method gives very similar results. For comparison, we use Matlab's 
ode45 method, that is explicit, fourth-order, and uses step size control. The 
ode45 method is not volume preserving unless the solution is computed to ma- 
chine accuracy. The standard implementation of ode45 (with stepsize control) 
becomes unstable at T = 2.2887e-|-2 (see illustration bottom right). In the bot- 
tom left picture we display the integration up to T = 220. To obtain a result 
similar to the volume-preserving method, it is necessary to decrease the error 
tolerance. The top right plot illustrates the results by ode45 with options set 
to RelTol=le-6 (default value le-3). 

Figure [3] shows the histogram for the stepsize chosen by ode45 to meet the 
required tolerance. For the smallest error tolerance, the average step size is 
around h = 0.01, which justifies the choice for our explicit volume-preserving 
methods. Numerical tests revealed that our methods were stable up to a choice 
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Figure 3: Histogram of the step sizes used by ode45 to meet the prescribed tol- 
erance for the quadratic volume-preserving flow (20). Dark gray: RelTol=le-6. 
Light gray: RelTol=le-3. 



h ~ 0.05. We performed also long time integration, T — 100000, with h = 0.05, 
and the solution still stayed bounded for the volume-preserving method. 

As far as computational time is concerned, our methods have the advantage 
of being explicit. The speedup, with respect to ode45 implemented with the 
option RelTol=le-6 for stability, is approximatly 3.2 for the first-order volume- 
preserving method and 3.15 for the second-order method, over an average of 
100 runs with the same initial condition as above, see Table [l] 





ode45 


vol. pres., order 1 


vol. pres., order 2 


cpu time [sec] 


3.8687 


1.2099 


1.2260 


speed up 


1 


3.1976 


3.1557 



Table 1: CPU time (in seconds) for ode45 with RelTol=le-6 and the explicit 
volume-preserving splitting methods. The values, for indication only, are com- 
puted as an average of 100 runs on an standard laptop (MacBook Pro) with 
initial condition Xq = [0, 0, 0.96]"^ and T = 500. 



6.2 Cubic Stokes flow 

We consider a cubic volume-preserving flow which has been used for the in- 
vestigation of the kinematics of bounded steady Stokes flows, in particular, to 
investigate the streamlines inside a neutrally buoyant spherical drop immersed 
in a genera linear flow (14] . This system is described by the divergence- free 
differential equation 

i = ^[{br'^ - 3)£:x - 2x(x^£;x)] + x 

where — x^x, w is the vorticity vector and E is the symmetric traceless 
rate-of-strain tensor of the external motion, that, in dimensionless terms, takes 
the form 

/ l/(l + a) 

E = i a/{l + a) 

V 0-1 
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Figure 4: The cubic Stokes flow (21) with w = (integrable case). Simulation 
of trajectories for T — 20000 for various initial conditions (the circles on the 
line denote the starting value). Top left: volume-preserving method, order 2, 
implemented with stepsize h — 0.01. Top right: volume-preserving method, 
order 1, implemented with stepsize h — 0.01. Bottom left: ode45, using step 
size control, with RelTol=le-6. Bottom right: ode45, default implementation 
(RelTol=le-3). 



where a = E22/E11, and x is the usual cross product. 

In Cartesian coordinates, the dynamical system for investigating fluid par- 
ticle motion takes the dimensionless form: 



Xi = 

±2 = 

±3 = 



(5r2-3)ffJ 
-(5r2- 3)2:3 



2^1 (iSs 

2X2 

2a;3 I i+„ 



l+a 
l+a 
l+a 



+ ^{W2X3 - W3X2) 
+ ^{W3X1 - W1X3) 
+ 5(^13^2 " W2X1). 



(21) 



The case w = is integrable. The generic case v^r 7^ is not integrable, but also 
in this case, solutions with initial condition in the unit sphere are bounded to 
the unit sphere, see [14; for more details. 



In the following experiments, the cubic Stokes flow (21) was computed nu- 



merically by volume-preserving methods and Matlab's ode45 method. Figure 
[4] shows the comparison of the trajectories by the volume-preserving methods 
(order 1 and 2) and Matlab's ode45 method with options RelTol=le-6 and 
RelTol=le-3 for w = 0. The volume-preserving methods preserve the qual- 
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itative behaviour for various initial conditions (denoted by small circles) in 
Figure |4] However, Matlab's ode45 method, with standard options control, 
has a dissipative behaviour. To obtain the same visual result as the volume- 
preserving methods, one has to reduce the tolerance on the error, for instance 
set RelTol=le-6. See the picture at the bottom of Figure |4] for more details. 

In the numerical investigation presented below (see Figure [s]) , we study the 
effect of changing both the orientation and the magnititude of w. Without 
loss of generality, we will take w = {wi, 0,103). The orientation of w, mea- 
sured from the a^a-axis in the {xi,X3) plane, will be denoted by 0. To describe 
the three-dimensional particle paths, we present Poincare sections through the 
(xi, X3) plane. We compute the particle paths by solving (21 ) numerically using 
a volume-preserving method (order 2), Matlab's ode45 method with options 
RelTol=le-6 and standard setup (RelTol=le-3) respectively. For more detail 
about the experiments design, see [14] for reference. 

To start with, we examine the case — 0.2757r and \\w\\ ~ 1.5 illustrated in 
the first row of Figure [5] The left column is the Poincare section obtained by 
the volume-preserving method (order 2), while the middle and right columns are 
obtained by Matlab's ode45 method with options RelTol=le-6 and standard 
implementation (RelTol=le-3) respectively. The particle parth Poincare sec- 
tion, computed by our volume-preserving method and Matlab's ode45 method 
with options RelTol=le-6 are visually identical to the figure showed in [T3] (see 
Figure 3 in [l3] for detail) , while Matlab's ode45 method with standard options 
RelTol=le-3 gives a totally wrong section. 

We next study the Poincare section for a fixed orientation of the vorticity 
vector, but increasing magnitude of the vorticity. The second and third rows of 
Figure [5] show simulations for an orientation 0.27r from the xa-axis by increasing 
||w|| from 2.5 to 4.0. As the magnitude of w increases, islands are created and 
destroyed until ||w|| = 1.4. For larger magnitudes of the vorticity vector, the 
streamlines structure is again modified and becomes much more complex [14) . 

The last two rows in the Figure [5] show the cases where the vorticity orien- 
tation increases from 0.027r to OAtt with ||w|| ~ 2.0 fixed. 

For all the experiments illustrated in Figure [5] the simulations have been 
performed with T = 20000 and initial condition Xq = [-0.1689,0,-0.0437]'^. 
The volume-preserving method (order 2) is implemented with constant step 
size h = 0.01 (the last row with step size h = 0.005). Each row corresponds to 
different value of w as explained above. 

From Figure [5j we can immediately observe that the particle parths in the 
Poincare sections by the volume-preserving method behave qualitatively sim- 
ilarly to the exact solution, while, for non volume-preserving methods (like 
Matlab's ode45), stronger restrictions on the step size are required to obtain 
the desired dynamic result. 
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Figure 5: The cubic Stokes flow (21) with w ^ (non-integrable case). 



Simulation of trajectories for T = 20000 with initial condition xq = 
[-0.1689,0,-0.0437]^. Left column: volume -preserving method, order 2, step- 
size h = 0.01 (with exception of last row, for which h = 0.005). Middle col- 
umn: ode45 with option RelTol=le-6. Right column: ode45, standard imple- 
mentation (RelTol=le-3). Each row corresponds to a different value of w = 
(w^, 0, Wy). From top to bottom: ||w|| = 1.5, O = 0.2757r; ||w|| = 2.5, O = 0.27r; 
||w|| = 4.0, = 0.27r; |lw|| = 2.0, O = 0.02n; ||w|| = 2.0, 6> = 0.47r. 
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Figure 6: Laurent polynomials problem (22) with different initial values. Left: 
volume-preserving method, order 2, implemented with stepsize h = 0.001. Mid- 
dle: od.e45, using stepsize control, with options RelTol=le-8, RelTol=le-8, 
RelTol=le-10 and RelTol=le-6. Right: ode45, using step size control, with 
standard options RelTol=le-3. 



6.3 Laurent polynomial example 

In order to illustrate our method for Laurent-type volume-preserving fields, we 
consider the following example: 

12 1 2 > ^22) 
According to the approach presented in this paper, we split the above equation 



(22) into two vector fields Fi and F2, corresponding to the multi- index ji 



(-3,2) andj2 = (2,-3), 



X2 ^~ ^X-^ X2 X2, — ^X ^X'2 



We can integrate the vector fields Fi and F2 explicitly according to Theo- 
rem |3] By choosing different initial values around [—0.5689,0.0437], we com- 
pare the results obtained by our method to Matlab's ode45 method with options 
RelTol=le-8, RelTol=le-8, RelTol=le-10 and RelTol=le-6 and standard im- 
plentation RelTol=le-3 respectively. In Figure[6] we show the numerical results 
for the volume-preserving method and Matlab's ode45 method, for initial val- 
ues [-0.5689 + kA, 0.0437 + lAf, k,l = 0,1, A = 0.02 (a small square). The 
standard Matlab's ode45 implementation diverged for early values of T (right 
column). Thus, to obtain the correct dynamical behaviour, stronger restric- 
tions on the error had to be imposed (middle column): a relative tolerance 
RelTol=le-8 gave results similar to the volume-preserving method for two of 
the initial conditions, for another one, the error control needed to be sharpened 
(RelTol=le-10), while for the last one, it was sufficient to use RelTol=le-6. 

Figure[7]shows the trajectories obtained with initial values [0.0437, —0.5489]-^ 
and [-0.5489, 0.437]^ (opposite octants) using the volume-preserving method 
(order 2) with stepsize h = 0.001 and Matlab's ode45 methods with options 
RelTol=le-8 and standard implementation. Matlab's ode45 method with stan- 
dard implementation was again unstable, and the trajectories diverged before 
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Figure 7: Laurent polynomials problem ( 22 1 with initial conditions 
(0.0437,-0.5489), (-0.5489,0.437) in the different Octant in the plane. Left: 
volume-preserving method, order 2, implemented with stepsize h = 0.001. Mid- 
dle: ode45, using stepsize control, with options RelTol=le-8. Right: ode45, 
using stepsize control, with standard implementation RelTol=le-3. 



I RGLTol=1e-8 
|R9LTol=1e-3 



Figure 8: Histogram of the stepsizes used by ode45 to meet the prescribed 



tolerance for (22) in [0,10]. Gray: RelTol=le-8. Black: RelTol=le-3 



the end of the integration interval was reached. The trajectories of the system 



( 22 ) are very sensitive with respect to the initial condition (ill conditioned prob- 
lem) and small changes of initial values may cause tremendous difference in the 
actual trajectories. This type of problems is well known to affect the numerical 
methods by imposing constraints on the step size (small step size must be used 
to remain sufficiently close to the true trajectory). Also our volume-preserving 
methods were affected by this, as well as Matlab's ode45, but they required 
milder restrictions. 

Figure[8]shows the histogram for the stepsize chosen by ode45 with RelTol=le- 
and standard implementation (RelTol=le-3) respectively, to meet the required 
tolerance in [0,10]. For the case ode45 with options RelTol=le-8, the aver- 
age stepsize is the interval (0,0.005). Looser restrictions gave diverging solu- 
tions. For comparison, the volume-preserving methods could use h = 0.25 and 
still remain bounded. However, for longer time integrations, also the volume- 
preserving methods had to reduce the step size to avoid divergence. 
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7 Conclusion and further remarks 



In this paper we have presented new, expUcit, volume-preserving splitting meth- 
ods for Laurent polynomial divergence-free vector fields. The methods rely on 
a decomposition of the divergence into a monomial basis. For each monomial 
basis element, a corresponding divergence- free elementary vector field is iden- 
tified and integrated exactly (hence the solution is volume preserving). These 
elementary fields are composed by a splitting method, giving rise to a overall 
explicit volume-preserving method. 

The method proposed are a significant contribution to the understanding 
of volume-preserving methods. To our knowledge, there is no explicit volume- 
preserving method that can deal with arbitrary polynomial vector fields of ar- 
bitrary degree, with exception of the case of linear and quadratic vector fields. 

As far as stability is concerned, our methods were stable for sufficiently small 
step size, a feature that is common to explicit methods. An upper bound on the 
step size can be obtained by Prop[2j Beware that this bound can be unrealistic 
for long time numerical integration and further step size restriction might occur, 
especially if the solution is unbounded, which is often the case of divergence-free 
vector fields and and high degrees of nonlinearity. These problems are harder 
to integrate and most numerical methods will require step size reduction for 
stability. 

To conclude, let us sketch out a possible way to attach the problem of volume- 
preserving integration for generic functions: let {(/)i(x)} a set of basis functions. 
Consider the decomposition 

p(x) = V-f(x)=^K0,(x) 

i 

in the given basis. Then, for each i, because of the independence of the basis 
functions, one must have pi = 0. Thus, as long as we are able to recover all 
the terms in the given vector field contributing to the coefficient p^, we will 
automatically have a volume-preserving split vector field. The difficulty is to 
find suitable basis functions, in the sense that the corresponding vector fields 
must be easy to integrate exactly. One such example is the monomial basis 
for polynomials, x-^, for which we have demonstrated the existence of exact 
formulas. 
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